#ifndef __particles
#define __particles
#include <iostream>
#include <vector>
#include "type.h"
#include "physicalconstants.h"
//#include "stdio.h"
///单个粒子类信息，可以是任何种类粒子
class __Ptcls{
    public:
        __Ptcls(__Vect2<double>& pos, __Vect3<double>& vel, int& wt,\
                double& ene, double charge):position(pos), old_pos(pos), \
                velocity(vel), old_vel(vel), weight(wt), energy(ene), charge(charge){
                id = total_id;
                total_id ++;
                total_num ++;
                lorentz_gamma = 1.0 / sqrt(1.0 - (velocity * velocity) / (c_light_speed\
                            * c_light_speed));
                __Vect3<double> zero(0.0, 0.0, 0.0);
                efield = zero;
                bfield = zero;
        }
        __Ptcls(const __Ptcls& other){
            this->species = other.species;
            this->position = other.position;
            this->velocity = other.velocity;
            this->old_vel = other.old_vel;
            this->old_pos = other.old_pos;
            this->charge = other.charge;
            this->id = other.id;
            //cout<<"in ptcls copy func "<<other.id<<endl;
            this->weight = other.weight;
            this->energy = other.energy;
            this->efield = other.efield;
            this->bfield = other.bfield;
            this->mass = other.mass;
            this->lorentz_gamma = other.lorentz_gamma;
            total_num ++;
        }
        __Ptcls(){
            id = total_id;
            total_id ++;
            total_num ++;
            __Vect3<double> zero(0.0, 0.0, 0.0);
            velocity = zero;
            old_vel = zero;
            lorentz_gamma = 1.0;
            energy = 0.0;
            efield = zero;
            bfield = zero;
        }

        __Ptcls(const Particles& myptcl);
    public:
        /// current max id 
        static long total_id;
        /// current box total num
        static long total_num;
        /// 粒子的标号
        long id;
        /// 粒子的种类：0电子，1：离子，2：正电子，3：光子
        int species;
        /// 粒子坐标,全局绝对
        __Vect2<double> position;
        /// 前一步的位置
        __Vect2<double> old_pos;
        /// 粒子速度
        __Vect3<double> velocity;
        /// velocity previous step
        __Vect3<double> old_vel;
        /// 粒子所在位置电场
        __Vect3<double> efield;
        /// 粒子所在位置磁场
        __Vect3<double> bfield;
        /// 粒子权重
        long weight;
        /// 粒子能量，没有乘以权重
        double energy;
        /// 粒子电荷，没有乘以权重
        double charge;
        //double charge;
        /// 粒子质量，没有乘以权重
        //double mass;
        double mass;
        double lorentz_gamma;
        //double lorentz_gamma;
        //double od_factor;
        double od_factor;
        //double rod_factor;
        double rod_factor;
        double qed_factor;
        //double qed_factor;
        ~__Ptcls(){
            total_num --;
        }
        friend ostream& operator<<(ostream& out, const __Ptcls& myptcl)
        {
            out<<myptcl.id<<"\t"<<myptcl.position<<"\t"<<myptcl.old_pos<<"\t"<<myptcl.velocity<<"\t"\
                <<myptcl.weight<<"\t"<<myptcl.energy<<"\t"<<myptcl.lorentz_gamma<<"\t"<<myptcl.charge<<"\t"\
                <<myptcl.mass;
            return out;
        }
        /// 粒子推动
        void advance(double& delta_t);
        /// boris转动进行电磁加速
        void boris(double& delta_t);
        /// 粒子根据速度推动一步
        void lorentz(double& delta_t);
        /// 设置粒子感受到的场
        void set_field(const __Vect3<double>& oefield, const __Vect3<double>& obfield){
            efield = oefield;
            bfield = obfield;
        }

};



void __Ptcls::advance(double& delta_t){
    boris(delta_t);
    lorentz(delta_t);
}

void __Ptcls::lorentz(double& delta_t){
    old_pos = position;
    __Vect2<double> vel(velocity + old_vel);
    position = position + vel * delta_t * 0.5;
#ifdef __Debug
    if(std::isnan(energy)) cout<<"nan energy in lorentz "<<endl;
#endif
    //getchar();
    old_vel = velocity;
}


void __Ptcls::boris(double& delta_t){
    //cout<<" pid "<<id<<" ps "<<species<<endl;
    if (mass < e_mass * 0.1)
    {
        return;
    }
    /***
    if(species == 1)
    {
        cout<<*this<<endl;
    }
    ***/
    __Vect3<double> zero(0, 0, 0);
    lorentz_gamma = 1.0 / sqrt(1.0 -(velocity * velocity) / c_light_speed / c_light_speed);
    double qtm2 = delta_t * charge * 0.5 / mass;
    __Vect3<double> uinit = lorentz_gamma * velocity;
    __Vect3<double> uminus = uinit + qtm2 * efield;
    
    lorentz_gamma = sqrt(1.0 + (uminus * uminus) / c_light_speed / c_light_speed);
    __Vect3<double> tvec = qtm2 * bfield / lorentz_gamma;
    __Vect3<double> uslash = uminus + Cross(uminus, tvec);
    __Vect3<double> svec = 2.0 * tvec / (1.0 + tvec * tvec);
    __Vect3<double> uplus = uminus + Cross(uslash, svec);
    
    __Vect3<double> ufinal = uplus + qtm2 * efield;
    lorentz_gamma = sqrt(1.0 + (ufinal * ufinal) / c_light_speed / c_light_speed);
    velocity = ufinal / lorentz_gamma;
    energy = (lorentz_gamma - 1.0) * mass * c_light_speed * c_light_speed;
    if(energy < 0)
    {
        cout<<" negative energy !"<<endl;
    }
    
    if(Mod(velocity)>c_light_speed)
    {
        cout<<" velocity error "<<efield<<endl;
    }

}



/***
 * how to cycle this ptclslist?
 *      1st declare an list and
 * typedef vector<vector<__Electron> > mylist;
 *      2nd declare an iterator
 * vector<__Electron>::iterator myiter;
 *      3rd use this method to cycle and
 * for (myiter = myelist.begin(); myiter != myelist.end(); myiter ++)
 * {
 *      cout<<"id of my ptcls"<<myiter->id<<endl;
 *      }
 *
 * I think i should write a iterator, so that it will be more flexiable
 *
 *
 *      ***/


 ///一个同类粒子的粒子束，隶属于一个网格，也可以不是
class __Ptclslist{
    public:
        /// 粒子类成员
        vector<__Ptcls> member;
        // use this func you can call by obj.iter() to get a iter
        // then use it to find any particle you want
        /// 粒子类的一个迭代器，一般不要调用，除非你懂
        typename vector<__Ptcls>::iterator iter;
        /// 当前粒子种类
        int species;
        /// 网格长宽
        __Vect2<double> mylength;
        // 0 e- 1 i+ 2 e+ 4 gamma
        __Ptclslist(){}
        __Ptclslist(const __Vect2<double>& mylength, int sp):mylength(mylength), species(sp){}
        __Ptclslist(const __Ptclslist& other){
            this->member = other.member;
            this->species = other.species;
            this->mylength = other.mylength;
        }
        /// 压缩空间
        void squeeze()
        {
            if (member.capacity() > 2 * member.size())
            {
                vector<__Ptcls>(member).swap(member);
            }
        }
        /// 利用此方法对当前的list进行压缩,如光子是需要压缩的.当光子的能量差异在delta_e之内时会发生压缩.当然可以指定空间距离
        void merge(const double& delta_e, const double& delta_v, const int& maxppc);
        // use this push func, different species ptcls cannot be
        // stored in same vector, since the are different __Ptcls(ty
        // pe)
        /// 将一个粒子压入当前list
        void push(const __Ptcls& other)
        {
            member.push_back(other);
            //cout<<"in pushing -----------------------------------"<<endl;
            //cout<<other.lorentz_gamma<<endl;
            //cout<<member.size()<<endl;
            //cout<<"pushin over-----------------------------"<<endl;
            //getchar();
        }
        // delete particle by offset, i.e. index, this is quicker
        // than other method
        // how to call?
        // just pop(n) n is the offset
        /// 删除offset处的粒子，一般不会主函数调用，不要使用，除非你懂
        void pop(int offset)
        {
            //cout<<"int plist offset"<<endl;
            if (offset < member.size())
            {
                iter = member.begin();
                std::advance(iter, offset);
                member.erase(iter);
            }
        }
        /// 清除当前list
        void clear()
        {
            clear_vec(member);
            //member.clear();
        }
        /// 返回当前list中粒子个数
        int size(){
            return member.size();
        }
        // use this func, you can get the begin of member, no inter
        // need to care about what kind of its type
        typename vector<__Ptcls>::iterator begin(){
            return member.begin();
        }
        typename vector<__Ptcls>::iterator end(){
            return member.end();
        }
        /// 返回offset处粒子信息
        __Ptcls& get_ptcl(int offset)
        {
            std::advance(iter, offset - 1);
            return *(iter);
        }
        /// 返回该list的能量
        double get_sp_ek(const int& sp);
        ~__Ptclslist(){
            clear_vec(member);
            //member.clear();
            this->squeeze();
        }
        /// 推动该list中所有粒子一步
        void advance(double delta_t, vector<__Ptcls>& lists, \
                const __Vect2<double>& xymin, const __Vect2<double>& xymax,\
                const int& myrank){
            if(member.size() == 0)
            {
                return;
            }
            int i = member.size() - 1;
            while(i >= 0){
                member[i].advance(delta_t);
                if((member[i].position.member[0] < xymin.member[0] || \
                            member[i].position.member[0] > xymax.member[0] || \
                            member[i].position.member[1] < xymin.member[1] || \
                            member[i].position.member[1] > xymax.member[1]))
                {
                    lists.push_back(member[i]);
                    pop(i);
                }
                i --;
            }
        }
        /// 输出每个粒子
        void showme(){
            for(int i = 0; i< member.size();i ++){
                cout<<member[i]<<endl;
            }
        }
        /// 通过重载可以直接复制了
        __Ptclslist operator = (const __Ptclslist& other){
            this->member = other.member;
            this->species = other.species;
        }
        /// 与cellinfo中相同
        double weightfunc(const __Vect2<double>& mycoords, const __Vect2<double>& other, \
                const __Vect2<double>& cell_length);
        double weightfunc(const __Vect2<double>& mycoords, const __Vect2<double>& other);
        double weight_to_coords_rho(const __Vect2<double>& coords, const __Vect2<double>&\
        cell_length);
        double weight_to_coords_j(const int& sp, const __Vect2<double>& coord, \
                const __Vect2<double>& cell_length);
        // 这里的sp指的是0---jx， 1--jy，2--jz
};

double __Ptclslist::weightfunc(const __Vect2<double>& mycoords, const __Vect2<double>& other, \
        const __Vect2<double>& cell_length){
    __Vect2<double> temp(0, 0);
    temp = (mycoords - other).Abs();
    //cout<<"temp "<<temp<<endl;
       if(cell_length > temp){
           return ((cell_length.member[0] - temp.member[0])\
                   * (cell_length.member[1] - temp.member[1]))\
               / (cell_length.member[0] * cell_length.member[1]);
       }
    else {
        return 0;
    }
}



double __Ptclslist::weight_to_coords_rho(const __Vect2<double>& coords, const __Vect2<double>&\
        cell_length){
    double rho = 0;
    if(member.size() != 0){
        for(int i = 0; i < member.size(); i ++){
            rho += member[i].weight * member[i].charge * \
                   weightfunc(member[i].position, coords, cell_length);

        }
    }
    return rho;
}

double __Ptclslist::weight_to_coords_j(const int& sp, const __Vect2<double>& coord, \
        const __Vect2<double>& cell_length){
    double rets = 0;
    if(member.size()!=0){
        for(int i = 0; i < member.size(); i ++){
            rets += member[i].weight * member[i].velocity.member[sp] *\
                    weightfunc(member[i].position, coord, cell_length) * member[i].charge;
        }
    }
    return rets;
}

double __Ptclslist::get_sp_ek(const int& sp)
{
    if(sp != species)
    {
        return 0;
    }
    double retek = 0;
    for(int i = 0; i < member.size(); i ++)
    {
        retek += member[i].weight * member[i].energy;
    }
    return retek;
}

void __Ptclslist::merge(const double& delta_e, const double& delta_v, const int& maxppc)
{
#ifdef __Debug
    ofstream myout;
    myout.open("dump.txt", ios::app);
#endif
    int nid = 0;
    __Vect3<double> dvv;
    //cout<<" merging "<<member.size()<<" -------> ";
    while(nid < member.size())
    {
        //cout<<member[nid].velocity<<"\t"<<member[nid].energy<<" is nid "<<member[nid].id<<endl;
        for(auto p = member.begin() + nid; p != member.end(); p ++)
        {
            //cout<<p->velocity<<"\t"<<p->energy<<" is p "<<p->id<<endl;
            dvv = member[nid].velocity - p->velocity;
            double delta_energy = abs(member[nid].energy - p->energy);
            double delta_velocity = Mod(dvv);
            /// 能量上加上1ev,防止粒子的能量为0, 速度上加1m/s,防止粒子完全静止
            if(delta_energy / (member[nid].energy + 1e-19) <= delta_e && \
                    (delta_velocity / (Mod(member[nid].velocity) + 1.0) < delta_v))
            //if(delta_energy / member[nid].energy <= delta_e && delta_energy / p->energy <= delta_e && \
                    (delta_velocity / Mod(member[nid].velocity) < delta_v && Mod(dvv) / \
                        Mod(p->velocity) < delta_v))
            {
                member[nid].energy = (member[nid].energy * member[nid].weight + p->energy * p->weight);
                member[nid].velocity = (member[nid].velocity * member[nid].weight + p->velocity * p->weight);
                member[nid].position = (member[nid].position + p->position) * 0.5;
                member[nid].weight += p->weight;
                member[nid].energy /= member[nid].weight;
                member[nid].velocity = member[nid].velocity / member[nid].weight;
                /// 删除被merge的粒子 
                p = member.erase(p);
                p --;
            }
            if(member.size() <= maxppc) break;
        }
#ifdef __Debug
        myout<<member[nid]<<endl;
#endif
        nid ++;
    }
#ifdef __Debug
    myout.close();
#endif
    //cout<<" "<<member.size()<<endl;
    /// 进行压缩空间
#ifdef __Debug
    for(auto p = member.begin(); p != member.end(); p ++)
    {
        if(std::isnan(p->energy)) cout<<"nan member in merge "<<endl;
    }
#endif
    squeeze();
}


__Ptcls::__Ptcls(const Particles& myptcl)
{
    id = total_id;
    total_id ++;
    mass = myptcl.mass;
    charge = myptcl.charge;
    energy = myptcl.energy;
    species = myptcl.species;
    velocity = myptcl.velocity;
    lorentz_gamma = 1.0 + myptcl.energy / (myptcl.mass * c_light_speed * c_light_speed);
    
}

#endif
